The Bragg condition is given by $2d_{hkl}\sin\theta = \lambda$. Since we know that $d_{hkl} = \frac{1}{|\mathbf{g_{hkl}}|}$, we have
$\theta = \sin^{-1}(\frac{\lambda |\mathbf{g_{hkl}}|}{2})$
Also, because the reciprocal basis vectors are mutually orthgonal, we have the lengths simply as
$|\mathbf{g_{hkl}}| = (h^2/a^2 + k^2/a^2 + l^2 / c^2)^{1/2}$
In [1]:
import math
a = 3.9039
c = 4.1348
def calculate_two_theta(hkl):
two_theta = 2 * math.asin(1.541 * (hkl[0] / a ** 2 + hkl[1] / a ** 2 + hkl[2] / c ** 2) ** 0.5 / 2)
return math.degrees(two_theta)
two_thetas = {}
for hkl in [(0, 1, 1), (0, 1, 0), (0, 0, 1)]:
two_theta = calculate_two_theta(hkl)
print "Peak position for %s is %.1f degrees" % (hkl, two_theta)
#Store this as we will use these later.
two_thetas[hkl] = two_theta
To calculate the structure factor, we need the positions of all atoms in the unit cell. Note that one of the oxygen has Wyckoff symbol 2c. From the International Tables, we find that the orbit for this site is given by (1/2, 0, z) and (0, 1/2, z). So there are 5 atoms in the cell with positions:
In [2]:
elements = ["Pb", "Ti", "O", "O", "O"]
coords = [
[0, 0, 0],
[0.5, 0.5, 0.5281],
[0.5, 0, 0.6130],
[0, 0.5, 0.6130],
[0.5, 0.5, 0.1339]
]
#We will also set up the parameters for the atomic scattering factors
asf_params = {
"Pb": [82, 3.510, 52.914, 4.552, 11.884, 3.154, 2.571, 1.359, 0.321],
"Ti": [22, 3.565, 81.982, 2.818, 19.049, 1.893, 3.590, 0.483, 0.386],
"O": [8, 0.455, 23.780, 0.917, 7.622, 0.472, 2.144, 0.138, 0.296]
}
The atomic scattering factor is given by:
$f(s) = Z - 41.78214 s^2 \sum(a_i e ^ {-b_i s^2})$
In [3]:
wavelength = 1.541 / 10 #Convert to nm
#Let's write a function to calculate the atomic scattering factor for a given element and angle.
def get_atomic_scattering_factor(el, two_theta):
params = asf_params[el]
theta = math.radians(two_theta / 2)
s = math.sin(theta) / wavelength
return params[0] - 41.78214 * s ** 2 * sum([params[2*i + 1] * math.exp(-params[2*i + 2] * s ** 2) for i in xrange(4)])
The structure factor is given by the sum of the atomic scattering factor times the position factor.
$F_{hkl} = \sum f_j(s) \exp(2\pi i (hx_j + ky_j + lz_j))$
and the intensity is simply given by the modulus squared of the structure factor.
In [4]:
import cmath # We need complex representation for this.
I = {}
for hkl in [(0, 1, 1), (0, 1, 0), (0, 0, 1)]:
F_hkl = 0
two_theta = two_thetas[hkl]
h, k, l = hkl
for el, coord in zip(elements, coords):
x, y, z = coord
F_hkl += get_atomic_scattering_factor(el, two_theta) * cmath.exp(2j * math.pi * (h * x + k * y + l * z))
print "Structure factor for %s = %s" % (hkl, str(F_hkl))
I_hkl = (F_hkl * F_hkl.conjugate()).real
I[hkl] = I_hkl
print "Intensity for %s = %.0f" % (hkl, I_hkl)
In [5]:
lp = {}
for hkl in [(0, 1, 1), (0, 0, 1)]:
theta = math.radians(two_thetas[hkl] / 2)
lp[hkl] = (1 + math.cos(2 * theta) ** 2) / (math.sin(theta) ** 2 * math.cos(theta))
print "Lorentz factor for %s = %.3f" % (hkl, lp[hkl])
In [6]:
print "Before correction, we have I_011/I_001 = %.4f" % (I[(0, 1, 1)]/I[(0, 0, 1)])
print "After correction, we have I_011/I_001 = %.4f" % (I[(0, 1, 1)] * lp[(0, 1, 1)] / (I[(0, 0, 1)]* lp[(0, 0, 1)]))
Clearly, the Lorentz polarization factor has a significant impact on relative peak intensities! In this case, the relative peak intensity has decreased significantly. If you look at the experiment XRD pattern below, clearly the relative intensity is a lot closer to 0.4143 than 0.9305. The relative intensities are still wrong because there are other factors which have not been taken into account.
In [7]:
from IPython.display import Image, display
display(Image(filename=('./PbTiO3.png')))
This is not part of the question, but just to show the actual computed pattern after taking all factors into account, which shows much better agreement with the experimental pattern.
In [8]:
from pymatgen import read_structure
from pymatgen.analysis.diffraction.xrd import XRDCalculator
%matplotlib inline
c = XRDCalculator()
c.get_xrd_plot(read_structure("PbTiO3.cif"))
Out[8]:
In [8]: